Two avenues of analysis:
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS",'reshape2',
'influence.ME', 'sjPlot', "effects", 'visreg', 'viridis')
ipak(packages)
## ggplot2 car lme4 gsheet MASS
## TRUE TRUE TRUE TRUE TRUE
## reshape2 influence.ME sjPlot effects visreg
## TRUE TRUE TRUE TRUE TRUE
## viridis
## TRUE
# read in data -- google sheet called "Final Bee Resp Data"
url = 'https://docs.google.com/spreadsheets/d/1wT-QxSJJElhhJcIXlg2hDpFKHbLyFGuqNNj2iYvA8Vo/edit?usp=sharing'
bdta <- data.frame(gsheet2tbl(url))
summary(bdta)
## BeeID order Treatment Mstarved
## Length:60 Min. :1.0 Length:60 Min. :0.0767
## Class :character 1st Qu.:1.0 Class :character 1st Qu.:0.1213
## Mode :character Median :1.5 Mode :character Median :0.1378
## Mean :1.5 Mean :0.1411
## 3rd Qu.:2.0 3rd Qu.:0.1670
## Max. :2.0 Max. :0.1909
## M2 MF Itspan S
## Min. :0.1062 Min. :0.0969 Min. :0.003860 Min. :4.010e-05
## 1st Qu.:0.1808 1st Qu.:0.1762 1st Qu.:0.004400 1st Qu.:5.590e-05
## Median :0.2191 Median :0.2107 Median :0.004620 Median :6.130e-05
## Mean :0.2245 Mean :0.2154 Mean :0.004594 Mean :6.141e-05
## 3rd Qu.:0.2646 3rd Qu.:0.2576 3rd Qu.:0.004870 3rd Qu.:6.750e-05
## Max. :0.3550 Max. :0.3440 Max. :0.005180 Max. :7.810e-05
## MetR freq amp L
## Min. : 5.507 Min. :163.9 Min. : 96.88 Min. :0.008691
## 1st Qu.: 8.935 1st Qu.:173.0 1st Qu.:114.16 1st Qu.:0.010570
## Median :10.666 Median :179.1 Median :125.96 Median :0.011017
## Mean :10.600 Mean :178.9 Mean :123.31 Mean :0.010911
## 3rd Qu.:12.297 3rd Qu.:185.3 3rd Qu.:133.61 3rd Qu.:0.011444
## Max. :16.475 Max. :194.7 Max. :144.78 Max. :0.012118
bdta <- within(bdta, {
MT = (M2 + MF)/2
load = MT - Mstarved
perLoad = load / Mstarved * 100
arcL = (0.75 * L) * (amp * (pi / 180))
U = arcL * freq * 2
freq2 = freq^2
arcL2 = arcL^2
frce = U^2 * S
# convert to factor variables
order = as.factor(as.character(order))
Treatment = as.factor(as.character(Treatment))
BeeID = as.factor(as.character(BeeID))
})
# convert from long to wide
newDF <- data.frame()
colsTocalc = c("order", "M2", "MF", "MetR", "freq", "amp", "load", "MT", "perLoad", "frce", "arcL2", "freq2", "U", "arcL")
for(varb in colsTocalc){
data_wide <- dcast(bdta, BeeID + Mstarved + S + Itspan + L ~
Treatment, value.var=c(varb))
colnames(data_wide)[6:7] = paste(varb, colnames(data_wide)[6:7], sep = "_")
if(varb == colsTocalc[1]){
newDF <- data_wide
}
else newDF <- merge(newDF, data_wide, all.y = TRUE)
}
head(newDF)
## BeeID Mstarved S Itspan L order_H order_L M2_H M2_L
## 1 E01 0.1577 6.46e-05 0.00464 0.01116843 2 1 0.2634 0.1895
## 2 E03 0.1732 6.92e-05 0.00487 0.01144364 1 2 0.3106 0.2105
## 3 E05 0.1755 7.06e-05 0.00497 0.01168730 2 1 0.2776 0.2006
## 4 E09 0.1135 4.88e-05 0.00440 0.00986884 1 2 0.2282 0.1425
## 5 E10 0.1269 5.91e-05 0.00462 0.01060036 1 2 0.2166 0.1460
## 6 E11 0.1350 5.88e-05 0.00464 0.01070018 1 2 0.2223 0.1599
## MF_H MF_L MetR_H MetR_L freq_H freq_L amp_H amp_L
## 1 0.2559 0.1812 12.54241 10.393483 186.9973 179.3308 128.7145 111.49875
## 2 0.3079 0.2023 11.78387 9.214622 174.9374 167.8412 144.7784 119.98953
## 3 0.2742 0.1975 12.20102 8.981412 178.8896 166.2385 125.6778 116.22053
## 4 0.2084 0.1390 9.47114 7.581104 191.8348 190.3535 139.3246 103.30703
## 5 0.2099 0.1422 11.03371 7.159840 186.1762 169.5451 128.8257 96.88003
## 6 0.2179 0.1579 10.23580 6.806696 183.0259 168.2223 128.7916 116.18262
## load_H load_L MT_H MT_L perLoad_H perLoad_L frce_H
## 1 0.10195 0.02765 0.25965 0.18535 64.64807 17.53329 0.003199482
## 2 0.13605 0.03320 0.30925 0.20640 78.55081 19.16859 0.003984230
## 3 0.10040 0.02355 0.27590 0.19905 57.20798 13.41880 0.003340857
## 4 0.10480 0.02725 0.21830 0.14075 92.33480 24.00881 0.002327019
## 5 0.08635 0.01720 0.21325 0.14410 68.04571 13.55398 0.002618299
## 6 0.08510 0.02390 0.22010 0.15890 63.03704 17.70370 0.002563877
## frce_L arcL2_H arcL2_L freq2_H freq2_L U_H
## 1 0.002208025 0.0003540922 0.0002657061 34967.99 32159.55 7.037583
## 2 0.002519159 0.0004703414 0.0003230668 30603.08 28170.67 7.587858
## 3 0.002467172 0.0003696773 0.0003161342 32001.50 27635.24 6.879020
## 4 0.001259710 0.0003239404 0.0001781022 36800.61 36234.44 6.905419
## 5 0.001228018 0.0003195386 0.0001807120 34661.58 28745.53 6.656039
## 6 0.001762569 0.0003254128 0.0002648146 33498.49 28298.73 6.603283
## U_L arcL_H arcL_L
## 1 5.846363 0.01881734 0.01630050
## 2 6.033576 0.02168736 0.01797406
## 3 5.911496 0.01922699 0.01778016
## 4 5.080722 0.01799834 0.01334549
## 5 4.558361 0.01787564 0.01344292
## 6 5.475004 0.01803920 0.01627313
newDF <- within(newDF, {
deltaPercLoad = perLoad_H - perLoad_L
avgPercLoad = (perLoad_H + perLoad_L) / 2
deltaMetR = MetR_H - MetR_L
deltaFrq2 = freq2_H - freq2_L
deltaArcL = arcL_H - arcL_L
deltaArcL2 = arcL2_H - arcL2_L
deltaFreq2Perc = deltaFrq2 / deltaPercLoad
deltaLoad = scale(load_H - load_L, center = TRUE, scale = FALSE)
dLoad_nonCent = load_H - load_L
deltaLoad2 <- deltaLoad^2
})
plot(newDF$deltaArcL2 ~ newDF$deltaFrq2)
# compare with Susie's calculations
b2 <- read.csv("~/Desktop/b2.csv")
b2 <- b2[!is.na(b2$deltaperload), ]
b2 <- b2[order(b2$BeeID, b2$Treatment), ]
plot(b2$AverageperLoad, newDF$avgPercLoad)
abline(0,1)
plot(b2$deltaload, newDF$deltaLoad)
abline(0,1)
# principle components
aa = prcomp(newDF[, c("Mstarved", "Itspan", "S", "L")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size
## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.9422 0.37234 0.22746 0.19332
## Proportion of Variance 0.9431 0.03466 0.01293 0.00934
## Cumulative Proportion 0.9431 0.97772 0.99066 1.00000
biplot(aa) # shows that all three size measurement are correlated
# note, I changed the signs of the predictions so that higher PC1 values
# correspond to bigger bees
p1 = -predict(aa)[,1]
# add PC1 scores to dataset
newDF$size_pc1 = p1
newDF$size_pc1_2 = newDF$size_pc1^2
# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(newDF[, c("Mstarved", "Itspan", "S", "L", "size_pc1")])
This is the model selection procedure that was used: 1. Fit a large model with all two-way interactions and squared terms 2. Remove non-significant predictors, starting with interactions, squared terms, and then main effects
# reformat order so that it is more interpretable
library(plyr)
newDF$order_1 <- mapvalues(newDF$order_H, from = c(2, 1), to = c("loadedSecond", "loadedFirst"))
# fit full model
m1 <- lm(deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 +
size_pc1_2 + deltaLoad2, data = newDF)
summary(m1)
##
## Call:
## lm(formula = deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 +
## size_pc1_2 + deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.326e-05 -1.951e-05 3.244e-06 1.635e-05 4.401e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.818e-05 1.300e-05 7.551 2.05e-07 ***
## deltaLoad 1.859e-03 1.178e-03 1.578 0.130
## size_pc1 -1.804e-05 1.263e-05 -1.428 0.168
## order_1loadedSecond 7.310e-07 1.978e-05 0.037 0.971
## size_pc1_2 -3.660e-06 4.905e-06 -0.746 0.464
## deltaLoad2 3.807e-02 4.366e-02 0.872 0.393
## deltaLoad:size_pc1 -2.300e-06 8.387e-04 -0.003 0.998
## deltaLoad:order_1loadedSecond -6.487e-04 2.225e-03 -0.292 0.774
## size_pc1:order_1loadedSecond 1.204e-05 2.099e-05 0.574 0.572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.03e-05 on 21 degrees of freedom
## Multiple R-squared: 0.575, Adjusted R-squared: 0.413
## F-statistic: 3.551 on 8 and 21 DF, p-value: 0.009387
m2 <- update(m1, .~. - deltaLoad:size_pc1)
anova(m1, m2)
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 +
## deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1 + size_pc1:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 1.9281e-08
## 2 22 1.9281e-08 -1 -6.9071e-15 0 0.9978
summary(m2)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + deltaLoad:order_1 + size_pc1:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.327e-05 -1.953e-05 3.263e-06 1.634e-05 4.402e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.817e-05 1.252e-05 7.840 8.24e-08 ***
## deltaLoad 1.860e-03 9.935e-04 1.872 0.0745 .
## size_pc1 -1.806e-05 7.886e-06 -2.290 0.0320 *
## order_1loadedSecond 7.471e-07 1.846e-05 0.040 0.9681
## size_pc1_2 -3.672e-06 1.745e-06 -2.104 0.0470 *
## deltaLoad2 3.800e-02 3.336e-02 1.139 0.2669
## deltaLoad:order_1loadedSecond -6.509e-04 2.022e-03 -0.322 0.7505
## size_pc1:order_1loadedSecond 1.209e-05 1.187e-05 1.019 0.3194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.96e-05 on 22 degrees of freedom
## Multiple R-squared: 0.575, Adjusted R-squared: 0.4397
## F-statistic: 4.251 on 7 and 22 DF, p-value: 0.004168
m3 <- update(m2, .~. - deltaLoad:order_1)
anova(m2, m3)
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1 + size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## size_pc1:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 22 1.9281e-08
## 2 23 1.9371e-08 -1 -9.0841e-11 0.1037 0.7505
summary(m3)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + size_pc1:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.414e-05 -1.856e-05 3.421e-06 1.715e-05 4.272e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.768e-05 1.218e-05 8.017 4.13e-08 ***
## deltaLoad 1.625e-03 6.585e-04 2.467 0.0215 *
## size_pc1 -1.728e-05 7.355e-06 -2.349 0.0278 *
## order_1loadedSecond 3.683e-06 1.573e-05 0.234 0.8170
## size_pc1_2 -3.816e-06 1.654e-06 -2.308 0.0304 *
## deltaLoad2 4.504e-02 2.467e-02 1.825 0.0809 .
## size_pc1:order_1loadedSecond 1.003e-05 9.807e-06 1.023 0.3169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.902e-05 on 23 degrees of freedom
## Multiple R-squared: 0.573, Adjusted R-squared: 0.4615
## F-statistic: 5.143 on 6 and 23 DF, p-value: 0.001756
m4 <- update(m3, .~. - size_pc1:order_1)
anova(m3, m4)
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 23 1.9371e-08
## 2 24 2.0253e-08 -1 -8.8175e-10 1.0469 0.3169
summary(m4)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.427e-05 -1.874e-05 -5.000e-09 1.783e-05 4.254e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.024e-04 1.128e-05 9.079 3.13e-09 ***
## deltaLoad 1.639e-03 6.590e-04 2.486 0.0203 *
## size_pc1 -1.180e-05 5.042e-06 -2.340 0.0279 *
## order_1loadedSecond -5.411e-07 1.519e-05 -0.036 0.9719
## size_pc1_2 -2.627e-06 1.178e-06 -2.231 0.0353 *
## deltaLoad2 2.613e-02 1.636e-02 1.597 0.1233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.905e-05 on 24 degrees of freedom
## Multiple R-squared: 0.5535, Adjusted R-squared: 0.4605
## F-statistic: 5.951 on 5 and 24 DF, p-value: 0.001027
m5 <- update(m4, .~. - deltaLoad2)
anova(m4,m5)
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 2.0253e-08
## 2 25 2.2406e-08 -1 -2.1529e-09 2.5511 0.1233
summary(m5)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.782e-05 -1.748e-05 9.200e-07 1.678e-05 4.278e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.074e-04 1.117e-05 9.623 6.93e-10 ***
## deltaLoad 2.095e-03 6.122e-04 3.422 0.00215 **
## size_pc1 -1.388e-05 5.018e-06 -2.767 0.01050 *
## order_1loadedSecond 2.783e-07 1.565e-05 0.018 0.98595
## size_pc1_2 -2.037e-06 1.152e-06 -1.768 0.08930 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.994e-05 on 25 degrees of freedom
## Multiple R-squared: 0.5061, Adjusted R-squared: 0.427
## F-statistic: 6.403 on 4 and 25 DF, p-value: 0.001086
m6 <- update(m5, .~. - size_pc1_2)
anova(m6, m5)
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 2.5207e-08
## 2 25 2.2406e-08 1 2.8009e-09 3.1251 0.0893 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m6)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.154e-05 -1.781e-05 -4.100e-07 1.849e-05 4.983e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.844e-05 1.034e-05 9.525 5.79e-10 ***
## deltaLoad 2.148e-03 6.359e-04 3.378 0.00231 **
## size_pc1 -1.321e-05 5.205e-06 -2.539 0.01745 *
## order_1loadedSecond 3.224e-06 1.618e-05 0.199 0.84364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.114e-05 on 26 degrees of freedom
## Multiple R-squared: 0.4443, Adjusted R-squared: 0.3802
## F-statistic: 6.929 on 3 and 26 DF, p-value: 0.001402
m7 <- update(m6, .~. - order_1)
anova(m7, m6)
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 2.5245e-08
## 2 26 2.5207e-08 1 3.8477e-11 0.0397 0.8436
summary(m7)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.002e-04 5.583e-06 17.940 < 2e-16 ***
## deltaLoad 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666
anova(m7, update(m7, .~. - deltaLoad)) # p-value for delta load
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ size_pc1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 2.5245e-08
## 2 28 4.5359e-08 -1 -2.0114e-08 21.512 8.055e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m7, update(m7, .~. - size_pc1)) # p-value for size
## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 2.5245e-08
## 2 28 3.4619e-08 -1 -9.3737e-09 10.025 0.003808 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m7)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.002e-04 5.583e-06 17.940 < 2e-16 ***
## deltaLoad 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666
# refit model with non-centered version of deltaLoad
# this model will have a different intercept, but same p-values for slopes
summary(update(m7, .~. - deltaLoad + dLoad_nonCent)) # final model for paper
##
## Call:
## lm(formula = deltaArcL2 ~ size_pc1 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.822e-05 3.460e-05 -1.683 0.10397
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## dLoad_nonCent 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666
par(mfrow = c(2,2))
plot(m7, which = 1:4) # no glaring violations, though there are a few fairly influential observations
par(mfrow = c(1,1))
car::vif(m7)
## deltaLoad size_pc1
## 1.841066 1.841066
summary(m7)
##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.002e-04 5.583e-06 17.940 < 2e-16 ***
## deltaLoad 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666
# calculate partial residuals for deltaLoad
# these are the residuals, minus the effect of detlaLoad
y <- residuals(m7, type = 'partial')[, "deltaLoad"]
# plot partial residuals with base R plotting
plot(x = newDF$deltaLoad, y = y)
# double check to make sure the slope for partial residual plots are the
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent))
##
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.584e-04 2.531e-05 -6.257 9.20e-07 ***
## newDF$dLoad_nonCent 2.059e-03 3.213e-04 6.409 6.14e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.003e-05 on 28 degrees of freedom
## Multiple R-squared: 0.5946, Adjusted R-squared: 0.5801
## F-statistic: 41.07 on 1 and 28 DF, p-value: 6.137e-07
# this is what the raw data look like
plot(x = newDF$deltaLoad, y = newDF$deltaFrq2)
# this package plots the partial residuals
crPlot(m7, variable = "deltaLoad")
# plot with ggplot2
theme_set(theme_classic())
# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaArcL, color = deltaLoad)) +
geom_point()
ggplot(newDF, aes(x= deltaLoad, y = deltaArcL, color = size_pc1)) +
geom_point()
# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) +
geom_point() +
labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on arcL^2 \n while subtracing affect of bee size") +
stat_smooth(method = 'lm', se = FALSE)
# fit full model
m1 <- lm(deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 +
size_pc1_2 + deltaLoad2, data = newDF)
summary(m1)
##
## Call:
## lm(formula = deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 +
## size_pc1_2 + deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3981.4 -1255.9 -170.3 1165.6 3552.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4195.97 909.39 4.614 0.00015 ***
## deltaLoad -106368.80 82390.69 -1.291 0.21073
## size_pc1 666.06 883.26 0.754 0.45917
## order_1loadedSecond -4376.19 1383.48 -3.163 0.00469 **
## size_pc1_2 -120.74 343.06 -0.352 0.72839
## deltaLoad2 -469149.49 3053199.27 -0.154 0.87935
## deltaLoad:size_pc1 6902.45 58657.40 0.118 0.90744
## deltaLoad:order_1loadedSecond 62840.01 155625.35 0.404 0.69045
## size_pc1:order_1loadedSecond 54.39 1467.95 0.037 0.97079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2119 on 21 degrees of freedom
## Multiple R-squared: 0.5048, Adjusted R-squared: 0.3162
## F-statistic: 2.676 on 8 and 21 DF, p-value: 0.03372
car::vif(m1) # some serious multicollinearity
## deltaLoad size_pc1 order_1
## 13.203325 19.003708 3.182227
## size_pc1_2 deltaLoad2 deltaLoad:size_pc1
## 18.740891 10.117446 25.115087
## deltaLoad:order_1 size_pc1:order_1
## 11.384208 23.301468
m2 <- update(m1, .~. - size_pc1:order_1)
anova(m1, m2)
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 +
## deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:size_pc1 + deltaLoad:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 94310333
## 2 22 94316499 -1 -6165.2 0.0014 0.9708
summary(m2)
##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + deltaLoad:size_pc1 + deltaLoad:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3973.1 -1259.6 -176.2 1149.9 3556.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4205.4 852.7 4.932 6.21e-05 ***
## deltaLoad -108212.1 64166.1 -1.686 0.10584
## size_pc1 695.7 367.7 1.892 0.07176 .
## order_1loadedSecond -4386.5 1324.0 -3.313 0.00316 **
## size_pc1_2 -109.9 176.0 -0.625 0.53864
## deltaLoad2 -417980.4 2660489.9 -0.157 0.87659
## deltaLoad:size_pc1 5130.1 33170.6 0.155 0.87850
## deltaLoad:order_1loadedSecond 66238.1 122843.3 0.539 0.59516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2071 on 22 degrees of freedom
## Multiple R-squared: 0.5048, Adjusted R-squared: 0.3472
## F-statistic: 3.203 on 7 and 22 DF, p-value: 0.01701
m3 <- update(m2, .~. - deltaLoad:size_pc1)
anova(m2, m3)
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:size_pc1 + deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 22 94316499
## 2 23 94419044 -1 -102545 0.0239 0.8785
summary(m3)
##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + deltaLoad:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3910.3 -1244.0 -168.5 1173.4 3503.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4190.99 829.34 5.053 4.09e-05 ***
## deltaLoad -108478.92 62767.13 -1.728 0.09734 .
## size_pc1 702.89 356.93 1.969 0.06108 .
## order_1loadedSecond -4434.40 1259.70 -3.520 0.00184 **
## size_pc1_2 -88.67 107.58 -0.824 0.41829
## deltaLoad2 -215246.26 2265485.45 -0.095 0.92513
## deltaLoad:order_1loadedSecond 61643.10 116639.62 0.528 0.60222
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2026 on 23 degrees of freedom
## Multiple R-squared: 0.5042, Adjusted R-squared: 0.3749
## F-statistic: 3.899 on 6 and 23 DF, p-value: 0.00785
m4 <- update(m3, .~. - deltaLoad:order_1)
anova(m3, m4)
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 23 94419044
## 2 24 95565635 -1 -1146591 0.2793 0.6022
summary(m4)
##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3667.1 -1308.7 -52.9 1317.9 3482.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.330e+03 7.749e+02 5.587 9.48e-06 ***
## deltaLoad -8.589e+04 4.527e+04 -1.897 0.069889 .
## size_pc1 7.352e+02 3.463e+02 2.123 0.044282 *
## order_1loadedSecond -4.794e+03 1.044e+03 -4.594 0.000117 ***
## size_pc1_2 -5.195e+01 8.090e+01 -0.642 0.526837
## deltaLoad2 -1.250e+06 1.124e+06 -1.112 0.277131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1995 on 24 degrees of freedom
## Multiple R-squared: 0.4982, Adjusted R-squared: 0.3937
## F-statistic: 4.766 on 5 and 24 DF, p-value: 0.00364
m5 <- update(m4, .~. - size_pc1_2)
anova(m4,m5)
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 95565635
## 2 25 97207799 -1 -1642165 0.4124 0.5268
summary(m5)
##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3591.2 -1388.0 -154.6 1391.9 3336.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4166.0 723.2 5.760 5.29e-06 ***
## deltaLoad -80714.2 44020.7 -1.834 0.078654 .
## size_pc1 732.5 342.2 2.141 0.042251 *
## order_1loadedSecond -4719.5 1024.9 -4.605 0.000104 ***
## deltaLoad2 -1475818.2 1054446.2 -1.400 0.173916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1972 on 25 degrees of freedom
## Multiple R-squared: 0.4896, Adjusted R-squared: 0.4079
## F-statistic: 5.995 on 4 and 25 DF, p-value: 0.00159
m6 <- update(m5, .~. - deltaLoad2)
anova(m6, m5)
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 104824695
## 2 25 97207799 1 7616896 1.9589 0.1739
summary(m6)
##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3735.3 666.5 5.605 6.87e-06 ***
## deltaLoad -105594.7 41007.6 -2.575 0.016065 *
## size_pc1 861.4 335.6 2.566 0.016382 *
## order_1loadedSecond -4717.7 1043.6 -4.521 0.000119 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared: 0.4496, Adjusted R-squared: 0.3861
## F-statistic: 7.08 on 3 and 26 DF, p-value: 0.001244
m7 <- update(m6, .~. - size_pc1)
anova(m7, m6) # p-values for size
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 131380648
## 2 26 104824695 1 26555953 6.5868 0.01638 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m8 <- update(m6, .~. - deltaLoad)
anova(m6, m8) # p-value for deltaLoad
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ size_pc1 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 104824695
## 2 27 131557512 -1 -26732816 6.6306 0.01607 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m9 <- update(m6, .~. - order_1)
anova(m6, m9) # p-value for order
## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 104824695
## 2 27 187214161 -1 -82389465 20.435 0.0001192 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit final model with non-centered deltaLoad
m10 <- update(m6, .~. - deltaLoad + dLoad_nonCent)
summary(m10) # final model for paper
##
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11857.7 3586.6 3.306 0.002766 **
## size_pc1 861.4 335.6 2.566 0.016382 *
## order_1loadedSecond -4717.7 1043.6 -4.521 0.000119 ***
## dLoad_nonCent -105594.7 41007.6 -2.575 0.016065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared: 0.4496, Adjusted R-squared: 0.3861
## F-statistic: 7.08 on 3 and 26 DF, p-value: 0.001244
par(mfrow = c(2,2))
plot(m10, which = 1:4) # no glaring violations
par(mfrow = c(1,1))
car::vif(m10) # vif is a little high
## size_pc1 order_1 dLoad_nonCent
## 3.056456 2.017003 3.643394
# calculate partial residuals for deltaLoad
# these are the residuals, minus the effect of detlaLoad
y <- residuals(m10, type = 'partial')[, "dLoad_nonCent"]
# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)
# double check to make sure the slope for partial residual plots are the
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent))
##
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8122 1631 4.980 2.93e-05 ***
## newDF$dLoad_nonCent -105595 20702 -5.101 2.11e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1935 on 28 degrees of freedom
## Multiple R-squared: 0.4816, Adjusted R-squared: 0.4631
## F-statistic: 26.02 on 1 and 28 DF, p-value: 2.106e-05
# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaFrq2)
# this package plots the partial residuals
crPlot(m10, variable = "dLoad_nonCent")
# plot with ggplot2
# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaFrq2, color = deltaLoad, shape = order_1)) +
geom_point()
ggplot(newDF, aes(x= deltaLoad, y = deltaFrq2, color = size_pc1, shape = order_1)) +
geom_point()
# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) +
geom_point() +
labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on freq^2 \n while subtracing affect of bee size and order") +
stat_smooth(method = 'lm', se = FALSE)
# partial residuals for order
y <- residuals(m10, type = 'partial')[, "order_1"]
ggplot(newDF, aes(x= order_1, y = y)) +
geom_boxplot() +
labs(x = "order", y = "partial residuals for order \n i.e. order effect on freq^2 \n while subtracing affect of bee size and load") +
stat_smooth(method = 'lm', se = FALSE)
# holding bee size and delta load constant, if a bee was loaded second, (confusing, huh?) then it would have a much lower freq^2 than if it was loaded first
summary(m10)
##
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11857.7 3586.6 3.306 0.002766 **
## size_pc1 861.4 335.6 2.566 0.016382 *
## order_1loadedSecond -4717.7 1043.6 -4.521 0.000119 ***
## dLoad_nonCent -105594.7 41007.6 -2.575 0.016065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared: 0.4496, Adjusted R-squared: 0.3861
## F-statistic: 7.08 on 3 and 26 DF, p-value: 0.001244
# make some predictions
predDF <- data.frame(size_pc1 = 0, order_1 = c("loadedSecond", "loadedFirst"), dLoad_nonCent = mean(newDF$dLoad_nonCent))
predDF$pred_deltaF2 <- predict(m10, predDF)
predDF
## size_pc1 order_1 dLoad_nonCent pred_deltaF2
## 1 0 loadedSecond 0.07692 -982.314
## 2 0 loadedFirst 0.07692 3735.337
mm1 <- lm(deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2, data = newDF)
car::vif(mm1)
## deltaFrq2 deltaArcL2 deltaLoad size_pc1 deltaLoad2
## 1.060059 1.852143 3.844322 2.679591 1.495113
summary(mm1)
##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad +
## size_pc1 + deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.11198 -0.64060 -0.01055 0.42457 2.97948
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.941e+00 6.182e-01 3.140 0.00444 **
## deltaFrq2 3.543e-04 6.888e-05 5.144 2.89e-05 ***
## deltaArcL2 -2.573e+03 5.900e+03 -0.436 0.66659
## deltaLoad 2.368e+01 1.937e+01 1.223 0.23337
## size_pc1 -5.482e-02 1.445e-01 -0.379 0.70776
## deltaLoad2 4.440e+02 5.114e+02 0.868 0.39381
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9233 on 24 degrees of freedom
## Multiple R-squared: 0.5937, Adjusted R-squared: 0.5091
## F-statistic: 7.014 on 5 and 24 DF, p-value: 0.0003629
mm2 <- update(mm1, .~. - size_pc1)
anova(mm1, mm2)
## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 20.460
## 2 25 20.582 -1 -0.12268 0.1439 0.7078
summary(mm2)
##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad +
## deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.06884 -0.60009 -0.03648 0.36122 3.09536
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.821e+00 5.223e-01 3.487 0.00182 **
## deltaFrq2 3.571e-04 6.730e-05 5.306 1.69e-05 ***
## deltaArcL2 -1.543e+03 5.147e+03 -0.300 0.76682
## deltaLoad 1.788e+01 1.168e+01 1.530 0.13849
## deltaLoad2 4.891e+02 4.888e+02 1.001 0.32665
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9074 on 25 degrees of freedom
## Multiple R-squared: 0.5913, Adjusted R-squared: 0.5259
## F-statistic: 9.041 on 4 and 25 DF, p-value: 0.0001166
mm3 <- update(mm2, .~. - deltaArcL2)
anova(mm3, mm2)
## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 20.656
## 2 25 20.582 1 0.07399 0.0899 0.7668
summary(mm3)
##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05850 -0.58135 -0.00991 0.37733 3.10956
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.682e+00 2.327e-01 7.228 1.12e-07 ***
## deltaFrq2 3.560e-04 6.602e-05 5.393 1.19e-05 ***
## deltaLoad 1.667e+01 1.077e+01 1.548 0.134
## deltaLoad2 4.422e+02 4.550e+02 0.972 0.340
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8913 on 26 degrees of freedom
## Multiple R-squared: 0.5898, Adjusted R-squared: 0.5425
## F-statistic: 12.46 on 3 and 26 DF, p-value: 3.065e-05
mm4 <- update(mm3, .~. - deltaLoad2)
anova(mm3, mm4)
## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 20.656
## 2 27 21.407 -1 -0.75053 0.9447 0.34
summary(mm4)
##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.824e+00 1.808e-01 10.086 1.18e-10 ***
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## deltaLoad 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06
mm5 <- update(mm4, .~. - deltaLoad)
anova(mm5, mm4) # p-value for load
## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 25.352
## 2 27 21.407 1 3.9448 4.9754 0.03421 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mm6 <- update(mm4, .~. - deltaFrq2)
anova(mm4, mm6) # p-value for deltafrq2
## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad
## Model 2: deltaMetR ~ deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 21.407
## 2 28 43.766 -1 -22.359 28.201 1.324e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# refit model with non-centered load
mm7 <- update(mm4, .~. - deltaLoad + dLoad_nonCent)
summary(mm7) # final model for paper
##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.777e-01 7.507e-01 0.237 0.8146
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06
car::vif(mm7)
## deltaFrq2 dLoad_nonCent
## 1.014369 1.014369
par(mfrow = c(2,2))
plot(mm7, which = 1:4) # no glaring violations
par(mfrow = c(1,1))
# calculate partial residuals for deltaLoad
# these are the residuals, minus the effect of detlaLoad
y <- residuals(mm7, type = 'partial')[, "dLoad_nonCent"]
# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)
# double check to make sure the slope for partial residual plots are the
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent))
##
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.6463 0.7371 -2.233 0.0337 *
## newDF$dLoad_nonCent 21.4030 9.3554 2.288 0.0299 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8744 on 28 degrees of freedom
## Multiple R-squared: 0.1575, Adjusted R-squared: 0.1274
## F-statistic: 5.234 on 1 and 28 DF, p-value: 0.02991
summary(mm7)
##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.777e-01 7.507e-01 0.237 0.8146
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06
# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaMetR)
# this package plots the partial residuals
crPlot(mm7, variable = "dLoad_nonCent")
crPlot(mm7, variable = "deltaFrq2")
summary(mm7) # final model for paper
##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.777e-01 7.507e-01 0.237 0.8146
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06
Refref: Do regression for Average percent loading vs. delta metabolic rate / 1% load etc. It’s the last page of the document Susie sent (3 different response variables).
Covariates: avg % load, order, bee size.
newDF <- within(newDF, {dmr_apl = deltaMetR / avgPercLoad})
mod1 <- lm(dmr_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1)
##
## Call:
## lm(formula = dmr_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.029235 -0.014142 -0.003264 0.008463 0.053034
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1404689 0.0176079 7.978 1.86e-08 ***
## avgPercLoad -0.0015015 0.0003086 -4.865 4.80e-05 ***
## order_1loadedSecond -0.0212737 0.0073765 -2.884 0.00779 **
## size_pc1 0.0012171 0.0019858 0.613 0.54527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01965 on 26 degrees of freedom
## Multiple R-squared: 0.6107, Adjusted R-squared: 0.5657
## F-statistic: 13.59 on 3 and 26 DF, p-value: 1.58e-05
mod2 <- update(mod1, .~. - size_pc1)
anova(mod1, mod2) # remove size_pc1
## Analysis of Variance Table
##
## Model 1: dmr_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: dmr_apl ~ avgPercLoad + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 0.010044
## 2 27 0.010189 -1 -0.00014511 0.3756 0.5453
summary(mod2) # final model for paper
##
## Call:
## lm(formula = dmr_apl ~ avgPercLoad + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030989 -0.013911 -0.001532 0.009837 0.050990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1431985 0.0168372 8.505 4.06e-09 ***
## avgPercLoad -0.0015573 0.0002914 -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551 0.0071702 -2.853 0.00822 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared: 0.605, Adjusted R-squared: 0.5758
## F-statistic: 20.68 on 2 and 27 DF, p-value: 3.578e-06
# get p-values for avgPercLoad and order
anova(mod2, update(mod2, .~. - order_1)) # p-value for order
## Analysis of Variance Table
##
## Model 1: dmr_apl ~ avgPercLoad + order_1
## Model 2: dmr_apl ~ avgPercLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 0.010189
## 2 28 0.013260 -1 -0.0030713 8.1384 0.008216 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod2, update(mod2, .~. - avgPercLoad)) # p-value for avg perc load
## Analysis of Variance Table
##
## Model 1: dmr_apl ~ avgPercLoad + order_1
## Model 2: dmr_apl ~ order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 0.010189
## 2 28 0.020965 -1 -0.010776 28.554 1.211e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# plot(I(deltaMetR) ~ avgPercLoad, data = newDF)
# plot(I(deltaMetR / avgPercLoad) ~ deltaMetR, data = newDF)
car::scatterplotMatrix(newDF[, c("avgPercLoad", "deltaMetR", "size_pc1")])
car::vif(mod1) # looks fine
## avgPercLoad order_1 size_pc1
## 1.114454 1.051699 1.116654
car::vif(mod2)
## avgPercLoad order_1
## 1.017216 1.017216
par(mfrow = c(2,2))
plot(mod2) # possible nonlinear trend in residuals
plot(mod1, which = 4)
# update model to add a non-linear term
newDF <- within(newDF, {avgPercLoad_cent = as.numeric(scale(newDF$avgPercLoad, center = TRUE, scale = FALSE))})
newDF$apl_cent2 <- with(newDF, avgPercLoad_cent^2)
m11 <- lm(dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m11)
##
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1 +
## size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030273 -0.010355 -0.003798 0.006817 0.055061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.247e-02 5.796e-03 9.053 2.29e-09 ***
## avgPercLoad_cent -1.553e-03 3.088e-04 -5.030 3.45e-05 ***
## apl_cent2 2.772e-05 2.280e-05 1.216 0.23544
## order_1loadedSecond -2.364e-02 7.564e-03 -3.125 0.00446 **
## size_pc1 1.321e-03 1.970e-03 0.671 0.50864
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01948 on 25 degrees of freedom
## Multiple R-squared: 0.6324, Adjusted R-squared: 0.5736
## F-statistic: 10.75 on 4 and 25 DF, p-value: 3.289e-05
car::vif(m11) # vif is much better with centered variables
## avgPercLoad_cent apl_cent2 order_1 size_pc1
## 1.135857 1.108436 1.126036 1.118753
m11a <- update(m11, .~. - size_pc1)
anova(m11, m11a) # remove size_pc1
## Analysis of Variance Table
##
## Model 1: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 25 0.0094835
## 2 26 0.0096541 -1 -0.00017058 0.4497 0.5086
m11b <- update(m11a, .~. - apl_cent2)
anova(m11a, m11b) # drop squared term
## Analysis of Variance Table
##
## Model 1: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1
## Model 2: dmr_apl ~ avgPercLoad_cent + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 0.0096541
## 2 27 0.0101893 -1 -0.00053522 1.4414 0.2407
summary(m11b)
##
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030989 -0.013911 -0.001532 0.009837 0.050990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0549528 0.0052157 10.536 4.58e-11 ***
## avgPercLoad_cent -0.0015573 0.0002914 -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551 0.0071702 -2.853 0.00822 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared: 0.605, Adjusted R-squared: 0.5758
## F-statistic: 20.68 on 2 and 27 DF, p-value: 3.578e-06
par(mfrow = c(2,2))
plot(m11b) # residuals look much better
plot(m11b, which = 4)
summary(m11b)
##
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030989 -0.013911 -0.001532 0.009837 0.050990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0549528 0.0052157 10.536 4.58e-11 ***
## avgPercLoad_cent -0.0015573 0.0002914 -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551 0.0071702 -2.853 0.00822 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared: 0.605, Adjusted R-squared: 0.5758
## F-statistic: 20.68 on 2 and 27 DF, p-value: 3.578e-06
## visualize model for deltametRate/avgPercLoading
par(mfrow = c(1,1))
plot(dmr_apl ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)
ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200),
order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))
ndf$apl_cent2 <- ndf$avgPercLoad_cent^2
ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(mod2, newdata = ndf), ndf)
ggplot(newDF, aes(x = avgPercLoad, y = dmr_apl)) +
geom_point(aes(color = order_1), , shape = 17) +
geom_line(data = predDF, aes(x = avgPercLoad, y = preds, color = order_1)) +
labs(x = "Average Load (% bodymass)",
y = "Change in Metabolic Rate (mL CO2 / hr) / \n Average Load (% bodymass)") +
scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) +
theme(legend.position = c(0.8,0.8))
ggsave("dmr_apl.pdf", width = 5, height = 4)
newDF <- within(newDF, {df2_apl = deltaFrq2 / avgPercLoad})
par(mfrow = c(1,1))
plot(df2_apl ~ avgPercLoad, newDF)
mod1_f <- lm(df2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_f)
##
## Call:
## lm(formula = df2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56.153 -25.697 3.939 22.568 51.988
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 221.4973 29.3034 7.559 5.04e-08 ***
## avgPercLoad -3.0335 0.5136 -5.906 3.14e-06 ***
## order_1loadedSecond -37.1107 12.2762 -3.023 0.00557 **
## size_pc1 -2.6032 3.3048 -0.788 0.43800
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.71 on 26 degrees of freedom
## Multiple R-squared: 0.667, Adjusted R-squared: 0.6286
## F-statistic: 17.36 on 3 and 26 DF, p-value: 2.15e-06
mod2_f <- update(mod1_f, .~. - size_pc1)
anova(mod1_f, mod2_f) # remove size_pc1
## Analysis of Variance Table
##
## Model 1: df2_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: df2_apl ~ avgPercLoad + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 27819
## 2 27 28482 -1 -663.86 0.6205 0.438
summary(mod2_f) # not final model for paper
##
## Call:
## lm(formula = df2_apl ~ avgPercLoad + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.479 -27.491 1.816 21.916 49.750
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 215.6591 28.1506 7.661 3.07e-08 ***
## avgPercLoad -2.9140 0.4873 -5.980 2.23e-06 ***
## order_1loadedSecond -38.8616 11.9881 -3.242 0.00315 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.48 on 27 degrees of freedom
## Multiple R-squared: 0.6591, Adjusted R-squared: 0.6338
## F-statistic: 26.1 on 2 and 27 DF, p-value: 4.904e-07
# get p-values for avgPercLoad
anova(mod2_f, update(mod2_f, .~. - order_1)) # p-value for order
## Analysis of Variance Table
##
## Model 1: df2_apl ~ avgPercLoad + order_1
## Model 2: df2_apl ~ avgPercLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 28482
## 2 28 39568 -1 -11086 10.508 0.003152 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(2,2))
plot(mod2_f) # non-linearity in residuals
plot(mod2_f, which = 4)
# update model to add a non-linear term
m22 <- lm(df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m22)
##
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1 +
## size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -35.430 -16.698 -1.678 19.332 43.880
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.70197 7.30757 4.886 5.01e-05 ***
## avgPercLoad_cent -3.27920 0.38928 -8.424 9.05e-09 ***
## apl_cent2 0.13215 0.02874 4.597 0.000106 ***
## order_1loadedSecond -48.37471 9.53594 -5.073 3.09e-05 ***
## size_pc1 -2.10870 2.48324 -0.849 0.403845
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.56 on 25 degrees of freedom
## Multiple R-squared: 0.8196, Adjusted R-squared: 0.7907
## F-statistic: 28.39 on 4 and 25 DF, p-value: 5.684e-09
car::vif(m22) # vif is much better with centered variables
## avgPercLoad_cent apl_cent2 order_1 size_pc1
## 1.135857 1.108436 1.126036 1.118753
m22a <- update(m22, .~. - size_pc1)
anova(m22, m22a) # remove size_pc1
## Analysis of Variance Table
##
## Model 1: df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 25 15074
## 2 26 15509 -1 -434.8 0.7211 0.4038
summary(m22a)
##
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.409 -16.168 -2.024 16.653 47.040
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.34581 7.22905 5.028 3.12e-05 ***
## avgPercLoad_cent -3.18454 0.37097 -8.584 4.60e-09 ***
## apl_cent2 0.13321 0.02856 4.664 8.17e-05 ***
## order_1loadedSecond -49.88054 9.31922 -5.352 1.33e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.42 on 26 degrees of freedom
## Multiple R-squared: 0.8144, Adjusted R-squared: 0.793
## F-statistic: 38.02 on 3 and 26 DF, p-value: 1.183e-09
par(mfrow = c(2,2))
plot(m22a) # residuals look much better
plot(m22a, which = 4)
summary(m22a) # final model for paper (though we could also report the non-centered values)
##
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.409 -16.168 -2.024 16.653 47.040
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.34581 7.22905 5.028 3.12e-05 ***
## avgPercLoad_cent -3.18454 0.37097 -8.584 4.60e-09 ***
## apl_cent2 0.13321 0.02856 4.664 8.17e-05 ***
## order_1loadedSecond -49.88054 9.31922 -5.352 1.33e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.42 on 26 degrees of freedom
## Multiple R-squared: 0.8144, Adjusted R-squared: 0.793
## F-statistic: 38.02 on 3 and 26 DF, p-value: 1.183e-09
## visualize model for deltametRate/avgPercLoading
ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200),
order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))
ndf$apl_cent2 <- ndf$avgPercLoad_cent^2
ndf$apl <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(m22a, newdata = ndf, se = TRUE), ndf)
par(mfrow = c(1,1))
plot(df2_apl ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)
ggplot(newDF, aes(x = avgPercLoad, y = df2_apl)) +
geom_point(aes(color = order_1), , shape = 17) +
geom_line(data = predDF, aes(x = apl, y = preds.fit, color = order_1)) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) +
labs(x = "Average Load (% bodymass)",
y = "Change in wingbeat freq^2 (hz^2) / \n Average Load (% bodymass)") +
scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) +
theme(legend.position = c(0.8,0.8))
ggsave("df2_apl.pdf", width = 5, height = 4)
newDF <- within(newDF, {da2_apl = deltaArcL2 / avgPercLoad})
par(mfrow = c(1,1))
plot(da2_apl ~ avgPercLoad, ylab = c("deltaArcL^2 / avgPercLoad"), newDF)
mod1_a <- lm(da2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_a)
##
## Call:
## lm(formula = da2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.023e-06 -5.865e-07 7.780e-08 4.151e-07 1.346e-06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.321e-06 6.130e-07 5.418 1.12e-05 ***
## avgPercLoad -2.018e-08 1.074e-08 -1.878 0.0717 .
## order_1loadedSecond -6.484e-07 2.568e-07 -2.525 0.0180 *
## size_pc1 1.764e-08 6.913e-08 0.255 0.8006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.843e-07 on 26 degrees of freedom
## Multiple R-squared: 0.3183, Adjusted R-squared: 0.2397
## F-statistic: 4.047 on 3 and 26 DF, p-value: 0.01737
mod2_a <- update(mod1_a, .~. - size_pc1)
anova(mod1_a, mod2_a) # remove size_pc1
## Analysis of Variance Table
##
## Model 1: da2_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: da2_apl ~ avgPercLoad + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 1.2174e-11
## 2 27 1.2204e-11 -1 -3.0501e-14 0.0651 0.8006
summary(mod2_a) # final model for paper
##
## Call:
## lm(formula = da2_apl ~ avgPercLoad + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.033e-06 -5.437e-07 9.477e-08 4.213e-07 1.332e-06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.361e-06 5.827e-07 5.768 3.91e-06 ***
## avgPercLoad -2.099e-08 1.009e-08 -2.081 0.0471 *
## order_1loadedSecond -6.365e-07 2.482e-07 -2.565 0.0162 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.723e-07 on 27 degrees of freedom
## Multiple R-squared: 0.3166, Adjusted R-squared: 0.266
## F-statistic: 6.255 on 2 and 27 DF, p-value: 0.005861
# get p-values for avgPercLoad and order
anova(mod2_a, update(mod2_a, .~. - order_1)) # p-value for order
## Analysis of Variance Table
##
## Model 1: da2_apl ~ avgPercLoad + order_1
## Model 2: da2_apl ~ avgPercLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 1.2204e-11
## 2 28 1.5178e-11 -1 -2.9742e-12 6.58 0.01619 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(mod2_a, update(mod2_a, .~. - avgPercLoad)) # p-value for avg perc load
## Analysis of Variance Table
##
## Model 1: da2_apl ~ avgPercLoad + order_1
## Model 2: da2_apl ~ order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 1.2204e-11
## 2 28 1.4161e-11 -1 -1.9567e-12 4.329 0.04708 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
par(mfrow = c(2,2))
plot(mod2_a)
plot(mod2_a, which = 4)
# visualize model
ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200),
order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))
ndf$apl_cent2 <- ndf$avgPercLoad_cent^2
ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(mod2_a, newdata = ndf, se = TRUE), ndf)
par(mfrow = c(1,1))
plot(da2_apl ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)
ggplot(newDF, aes(x = avgPercLoad, y = da2_apl)) +
geom_point(aes(color = order_1), shape = 17) +
geom_line(data = predDF, aes(x = avgPercLoad, y = preds.fit, color = order_1)) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) +
labs(x = "Average Load (% bodymass)",
y = "Change in arc length^2 (radians^2) / \n Average Load (% bodymass)") +
scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) +
theme(legend.position = c(0.8,0.8))
ggsave("da2_apl.pdf", width = 5, height = 4)
sessionInfo()
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plyr_1.8.4 viridis_0.3.4 visreg_2.3-0
## [4] effects_3.1-2 sjPlot_2.2.1 influence.ME_0.9-8
## [7] reshape2_1.4.2 MASS_7.3-45 gsheet_0.4.2
## [10] lme4_1.1-12 Matrix_1.2-8 car_2.1-4
## [13] ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyr_0.6.1 splines_3.3.2 merTools_0.3.0
## [4] modelr_0.1.0 shiny_1.0.0 assertthat_0.1
## [7] stats4_3.3.2 coin_1.1-3 yaml_2.1.14
## [10] backports_1.0.5 lattice_0.20-34 quantreg_5.29
## [13] digest_0.6.12 minqa_1.2.4 colorspace_1.3-2
## [16] sandwich_2.3-4 htmltools_0.3.5 httpuv_1.3.3
## [19] psych_1.6.12 broom_0.4.1 SparseM_1.74
## [22] haven_1.0.0 purrr_0.2.2 xtable_1.8-2
## [25] mvtnorm_1.0-5 scales_0.4.1 stringdist_0.9.4.4
## [28] MatrixModels_0.4-1 arm_1.9-3 tibble_1.2
## [31] mgcv_1.8-16 DT_0.2 TH.data_1.0-8
## [34] nnet_7.3-12 lazyeval_0.2.0 pbkrtest_0.4-6
## [37] mnormt_1.5-5 survival_2.40-1 magrittr_1.5
## [40] mime_0.5 evaluate_0.10 nlme_3.1-130
## [43] foreign_0.8-67 tools_3.3.2 multcomp_1.4-6
## [46] stringr_1.1.0 munsell_0.4.3 blme_1.0-4
## [49] grid_3.3.2 nloptr_1.0.4 htmlwidgets_0.8
## [52] labeling_0.3 rmarkdown_1.3 gtable_0.2.0
## [55] codetools_0.2-15 curl_2.3 abind_1.4-5
## [58] sjstats_0.7.1 DBI_0.5-1 sjmisc_2.2.1
## [61] R6_2.2.0 gridExtra_2.2.1 zoo_1.7-14
## [64] knitr_1.15.1 dplyr_0.5.0 rprojroot_1.2
## [67] readr_1.0.0 modeltools_0.2-21 stringi_1.1.2
## [70] parallel_3.3.2 Rcpp_0.12.9 coda_0.19-1
## [73] lmtest_0.9-34
Sys.time()
## [1] "2017-03-16 10:12:29 EDT"